#coding=utf-8
import math
import re
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from utils.mylib_utils import *

# 原位火箭：
# Malan_BRS_ATDW_14DS20220101_20240125140000_V01.01.DAT
def evaluate_uncertainty_rocket(L2_filepath, system_param, k=1, given_param="U"):
    print("评定参数1：", system_param)
    print("评定参数2：", given_param)

    L2_filename = os.path.basename(L2_filepath)
    info = L2_filename.split("_")
    # given_param = "T"               # 传参（测试）
    device_code = info[1]           # 设备编码：探空气球等
    data_type = info[2]
    data_L2 = parse_QQ_data(L2_filepath)    # 0高度,1气压,2温度,3密度,4风速,5风向
    # print(data_L2)
    if data_type == "ATDW":
        if given_param == "T":
            data_L2 = data_L2[:, (0, 2, 5)]
        elif given_param == "P":
            data_L2 = data_L2[:, (0, 1)]
        elif given_param == "D":
            data_L2 = data_L2[:, (0, 3)]
        elif given_param == "W":
            data_L2 = data_L2[:, (0, 4)]
        # elif given_param == "X":
            # data_L2 = data_L2[:, (0, 5)]

        param = given_param

    # print(data_L2)
    #return
    result_df = None
    # try:
    if param == "T":
        result_df = get_GUM_Rocket_T_uncertainty(
            data_L2,
            k,                                     # GUM参数k
            system_param["T_error"],               # 温度载荷的最大允许误差（输入），默认："(0,50]km:3.5;(50,70]km:0.35"
            system_param["k_square"],              # 系统所在海拔高度, km
            system_param["day_time"],              # 火箭释放时间：1白天，0晚上
            system_param["T_u2_delta_radio"],        # 辐射最大修正量（适用于白天day_time=1）
            system_param["delta_r_param"],         # 恢复系数误差（未用）
            system_param["Vz"],                    # 下降速度（按速度仿真结果读取，未提供）（T未用，P使用）
            system_param["delta_Vz"],               # 下降速度的精度(m/s)（T未用，P使用）
            system_param["cp"],                     # 空气的定压比热（单位：kJ/(kg·K)）（T未用）
            system_param["T_u3_delta_aerodynamic"],   # 气动最大修正值
            system_param["is_correct"],             # 是否已做修正
            system_param["u4_sluggish_param"],      # 迟滞影响量参数
            system_param["T_param_u5"],               # 探空仪传感结构引入的测量不确定度分量
        )
    elif param == "P":
        result_df = get_GUM_Rocket_P_uncertainty(
            data_L2,
            k,                                     # GUM参数k
            system_param["k_square"],              # 系统所在海拔高度, km
            system_param["P_error_param"],         # 压力传感器的最大允许误差
            system_param["P_u2_delta_T_features"], # 传感器温度特性半宽度
            system_param["P_param_u3"],            # 传感器气动效应引入的测量不确定度分量
            system_param["P_u4_p"],                # 空气密度（kg/m3）
            system_param["Vz"],                    # 下降速度（按速度仿真结果读取，未提供）（T未用，P使用）
            system_param["delta_Vz"],              # 下降速度的精度(m/s)（T未用，P使用）
        )
    elif param == "W":
        result_df = get_GUM_Rocket_W_uncertainty(
            data_L2,
            k,                                     # GUM参数k
            system_param["k_square"],              # 系统所在海拔高度, km
            system_param['delta_t'],               # 时间间隔测量的不确定度
            system_param['deltaV']                 # 风速系数
        )
    elif param == "D":
        print("原位探测火箭不评密度探测值..")
    # except:
    #     print('File Error: '+L2_filepath)
    return result_df


# 温度不确定度计算
def get_GUM_Rocket_T_uncertainty(
        data_L2,
        K,
        T_error = "(0,50]km:3.5;(50,70]km:0.35",      # 温度载荷的最大允许误差（输入）默认："(0,50]km:3.5;(50,70]km:0.35"
        k_square = 3,                                 # sqrt_k计算系数
        day_time = 1,                                 # 释放时间：1白天，0晚上
        u2_delta_radio = "(0,45]km:0.8;(45,55]km:1.0;(55,60]km:1.5;(60,70]km:2.3",    # 辐射最大修正量（适用于白天day_time=1）
        delta_r_param = "(0,40]km:4e-4;(40,50]km:1.6e-3;(50,60]km:8e-3;(60,70]km:10e-3",            # 恢复系数误差（未用）
        Vz = "8",                                      # 下降速度（按速度仿真结果读取，未提供）（未用）
        delta_Vz = "0.03",                            # 下降速度的精度(m/s)（未用）
        cp = 1.004,                                   # 空气的定压比热（单位：kJ/(kg·K)）（未用）
        u3_delta_aerodynamic = "(0,30]km:0.1;(30,40]km:0.6;(40,45]km:1.3;(45,50]km:3.0;(50,55]km:5.1;(55,60]km:6.3;(60,70]km:9.0",    # 气动最大修正值
        is_correct = 0,       # 是否已做修正，0：未修正，1：已修正
        # 方法一：
        # m = 0.1,            # 传感器敏感元件质量(g)
        # C = 0.14,           # 传感器热容量(kJ/（kg·K）)
        # A = 3.14*10-6,      # 传感器的表面积为0.5mm直径的球表面积，3.14×10-6m2；(m2)
        # 未完..
        # 方法二：
        u4_sluggish_param = "(0,50]km:0.1;(50,58]km:0.4;(58,60]km:0.9;(60,70]km:6.0",   # 迟滞影响量参数
        param_u5 = 0,       # 探空仪传感结构引入的测量不确定度分量
):
    # 获取高度列，温度列
    z = data_L2[:, 0] / 1000        # 高度单位转换为 km
    T = data_L2[:, 1]
    shape = np.shape(z)

    # 计算中使用的常量
    k_sqrt = math.sqrt(k_square)

    # 不确定度分量
    # 1.由温度传感器自身准确度引入的测量不确定度分量计算【完成】
    t_error_arr = parseParamStr(T_error)
    MPE = getParamByHeightV2(z, t_error_arr)
    u1 = np.full(shape, MPE / k_sqrt)
    u1 = u1.astype('float')

    # 2.由太阳辐射引入的测量不确定度分量计算
    if day_time == 1:
        delta_radio_param = parseParamStr(u2_delta_radio)
        delta_radio = getParamByHeightV2(z, delta_radio_param)
        u2 = np.full(shape, delta_radio / k_sqrt)
    else:
        u2 = 0
        u2 = np.full(shape, u2)
    u2 = u2.astype('float')

    # 3.由气动效应引入的测量不确定度分量计算
    delta_aerodynamic_range = parseParamStr(u3_delta_aerodynamic)
    delta_aerodynamic = getParamByHeightV2(z, delta_aerodynamic_range)
    u3 = np.full(shape, delta_aerodynamic / (2 * k_sqrt))
    u3 = u3.astype('float')

    # 4.由迟滞效应引入的测量不确定度分量计算
    # if is_correct == 0:
    sluggish_param = parseParamStr(u4_sluggish_param)
    sluggish = getParamByHeightV2(z, sluggish_param)
    u4 = np.full(shape, sluggish / (2 * k_sqrt))
    u4 = u4.astype('float')
    # else:
    # TODO：未考虑未修正的计算方法

    # 5.由探空仪传感结构引入的测量不确定度分量计算
    u5 = np.full(shape, param_u5)
    u5 = u5.astype('float')

    # print(        u1**2
    #               + u2**2
    #               + u3**2
    #               + u4**2
    #               + u5**2)

    # 计算合成标准不确定度
    # 合成标准不确定度
    u = np.sqrt(
        u1**2
        + u2**2
        + u3**2
        + u4**2
        + u5**2
    )
    # 扩展不确定度
    U = K * u
    # 不确定度
    U, u, u1, u2, u3, u4, u5 = np.round(
        (U, u, u1, u2, u3, u4, u5), 4
    )
    result = np.vstack((z, T, U, u, u1, u2, u3, u4, u5)).T
    header = [
        "高度(km)",
        "温度(K)",
        "U(K)_GUM扩展不确定度",
        "u(K)_标准不确定度",
        "u1(K)_由温度传感器自身准确度引入",
        "u2(K)_由太阳辐射引入",
        "u5(K)_由气动效应引入",
        "u7(K)_由迟滞效应引入",
        "u8(K)_由探空仪传感结构引入"
    ]
    result_df = pd.DataFrame(result, columns=header)
    return result_df


# 压力不确定度计算
def get_GUM_Rocket_P_uncertainty(
        data_L2,
        K,
        k_square = 3,               # sqrt_k计算系数
        P_error_param = "(1,1000]Pa:10;(1000,110000]Pa:50",      # 压力传感器的最大允许误差
        P_u2_delta_T_features = 0,      # 传感器温度特性半宽度
        P_param_u3 = 0,                 # 传感器气动效应引入的测量不确定度分量
        P_u4_p = 1.29,                  # 空气密度（kg/m3）
        Vz = "8",                       # 下降速度（按速度仿真结果读取，未提供）
        delta_Vz = "0.03",              # 下降速度的精度(m/s)
):
    # 获取高度列，压力
    z = data_L2[:, 0] / 1000        # 高度单位转换为 km
    P = data_L2[:, 1]
    shape = np.shape(z)

    # 计算中使用的常量
    k_sqrt = math.sqrt(k_square)

    delta_Vz = float(delta_Vz)
    Vz = float(Vz)

    # 不确定度分量
    # 1.由压力传感器自身准确度引入的测量不确定度分量计算
    P_error_range = parseParamStr(P_error_param)
    MPE = getParamByHeightV2(P, P_error_range)
    u1 = np.full(shape, MPE / k_sqrt)
    u1 = u1.astype('float')

    # 2.由传感器温度特性引入的测量不确定度分量计算
    u2 = np.full(shape, P_u2_delta_T_features / k_sqrt)
    u2 = u2.astype('float')

    # 3.由传感器气动效应引入的测量不确定度分量计算
    u3 = np.full(shape, P_param_u3)
    u3 = u3.astype('float')

    # 4.由风压引入的测量不确定度分量
    u4 = P_u4_p * (2 * delta_Vz * Vz + delta_Vz ** 2) / (2 * k_sqrt)
    u4 = np.full(shape, u4)
    u4 = u4.astype('float')

    # [计算合成标准不确定度]
    # 合成标准不确定度
    u = np.sqrt(
        u1**2
        + u2**2
        + u3**2
        + u4**2
    )
    # [扩展不确定度]
    U = K * u
    # 返回
    U, u, u1, u2, u3, u4 = np.round(
        (U, u, u1, u2, u3, u4), 4
    )
    result = np.vstack((z, P, U, u, u1, u2, u3, u4)).T
    header = [
        "高度(km)",
        "气压(Pa)",
        "U(K)_GUM扩展不确定度",
        "u(K)_标准不确定度",
        "u3(K)_由压力传感器自身准确度引入",
        "u4(K)_由温度特性引入",
        "u5(K)_由气动效应引入",
        "u14(K)_由风压引入"
    ]
    result_df = pd.DataFrame(result, columns=header)
    return result_df


# 风速不确定度计算
def get_GUM_Rocket_W_uncertainty(
        data_L2,
        K,
        k_square = 3,       # u1计算系数
        delta_t = 0.001,    # 时间间隔测量的不确定度
        deltaV = 0.3        # 风速系数
):
    # 获取高度列，温度列
    z = data_L2[:, 0] / 1000        # 高度单位转换为 km
    W = data_L2[:, 1]
    shape = np.shape(z)

    # 计算中使用的常量
    k_sqrt = math.sqrt(k_square)

    # 不确定度分量
    # 1.由火箭速度测量引入的测量不确定度分量计算
    u1 = delta_t * W

    # 2.由大气风速解算引入的测量不确定度分量计算（等同于气球风场不确定度分量u2的计算）
    deltaP = getSensorErrorValueByWindSpeed(W)
    u2 = deltaP / k_sqrt

    # 3.由气球风阻系数引入的测量不确定度分量计算
    u3_tmp = deltaV / k_sqrt
    u3 = np.full(shape, u3_tmp)

    # [计算合成标准不确定度]
    # 合成标准不确定度
    u = np.sqrt(
        u1**2
        + u2**2
        + u3**2
    )
    # [扩展不确定度]
    U = K * u
    # 精度处理
    U, u, u1, u2, u3 = np.round(
        (U, u, u1, u2, u3), 4
    )
    result = np.vstack((z, W, U, u, u1, u2, u3)).T
    header = [
        "高度(km)",
        "温度(K)",
        "U(K)_GUM扩展不确定度",
        "u(K)_标准不确定度",
        "u15(K)_由火箭速度测量引入",
        "u16(K)_由大气风速解算引入",
        "u17(K)_由火箭风阻系数引入",
    ]
    result_df = pd.DataFrame(result, columns=header)
    return result_df


# 根据风速判断参数
def getSensorErrorValueByWindSpeed(W):
    t_high_ratio_arr = [('-', 10, 1), (10, '-', 0.1)]
    ratio = np.full(np.shape(W), np.nan)
    for ratio_param in t_high_ratio_arr:
        if ratio_param[0] == '-':
            mask = (W < ratio_param[1])
            ratio[mask] = ratio_param[2]
        elif ratio_param[1] == '-':
            mask = (ratio_param[0] <= W)
            # ratio[mask] = ratio_param[2] * W
            ratio[mask] = ratio_param[2]
    return ratio


# 根据气压值判断压力传感器的最大误差
# 气压≥500hPa，△(P)≤2hPa
# 气压＜500hPa，△(P)≤1hPa
#
def getSensorErrorValueByPressure(P):
    t_high_ratio_arr = [('-', 500, 1), (500, '-', 2)]
    ratio = np.full(np.shape(P), np.nan)
    for ratio_param in t_high_ratio_arr:
        if ratio_param[0] == '-':
            mask = (P < ratio_param[1])
            ratio[mask] = ratio_param[2]
        elif ratio_param[1] == '-':
            mask = (ratio_param[0] <= P)
            ratio[mask] = ratio_param[2]
    return ratio


# 根据高度判断参数值
# 参数格式：
# t_high_ratio_arr = [(0, 11, -0.007), (11, 20, 0.00), (20, 32, 0.001)]
# t_high_ratio_arr可使用函数parseParamStr获取
def getParamByHeightV2(H, t_high_ratio_arr):
    ratio = np.full(np.shape(H), np.nan)
    for ratio_param in t_high_ratio_arr:
        mask = (ratio_param[0] <= H) & (H <= ratio_param[1])
        ratio[mask] = ratio_param[2]
    return ratio


# def getParamByHeightV3(z):
#     t_high_ratio = [(0, 11, -0.007), (11, 20, 0.00), (20, 32, 0.001)]
#
#     height_km = z / 1000
#     # height_km = np.array(z) / 1000
#     # height_km = z.apply(lambda x: x / 1000)
#     ratios = np.array(t_high_ratio)
#     # conditions = (height_km > ratios[:, 0]) & (height_km < ratios[:, 1])
#     # conditions = [((x > ratios[:, 0]) & (x < ratios[:, 1])).any() for x in zip(*height_km)]
#     conditions = np.any((height_km[:, None] > ratios[:, 0]) & (height_km[:, None] < ratios[:, 1]), axis=1)
#     # values = np.where(conditions, ratios[:, 2], 0)
#     # values = np.where(conditions, np.repeat(ratios[:, 2], len(height_km)), 0)
#     values = np.where(conditions, np.repeat(ratios[:, 2], ratios[:, 1] - ratios[:, 0]), 0)
#     return values.max()


# 解析分段函数参数为数组
# 输入："(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001"
# 输出：[(0, 11, -0.007), (11, 20, 0.00), (20, 32, 0.001)]
def parseParamStr(piecewise_fun_str):
    ret = []
    arr = re.split(";", piecewise_fun_str)
    for item in arr:
        arr1 = re.split(":", item)
        if len(arr1) >= 2:
            value = arr1[1]         # value = -0.007
            print(value)
            # 使用正则表达式进行匹配和提取
            flag = False
            pattern = r'\((.*?)\]'  # 匹配以括号开头，中间任意字符，以]结尾的部分。支持： (a, b]
            match = re.search(pattern, arr1[0])
            if match:
                flag = True
                condition_text = match.group(1)  # 获取匹配结果中的第一个捕获组（codition_text: 0,11; 11,20; 20,32）
                # print(condition_text)
            else:
                pattern = r'\((.*?)\)'  # 支持：(a, b)
                match = re.search(pattern, arr1[0])
                if match:
                    flag = True
                    condition_text = match.group(1)
            if flag:
                arr3 = re.split(",", condition_text)
                print(arr3)
                element = (float(arr3[0]), float(arr3[1]), float(value))
                ret.append(element)
            else:
                print("分段函数参数错误（level-2）！");
                continue
        else:
            print("分段函数错误（level-1）！")
            continue
    return ret


# 未使用（已废弃）
# 解析参数
# "(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001"
def getParamByHeight(height, param_str):
    height_km = float(height)/1000
    arr = re.split(";", param_str)
    for item in arr:
        arr1 = re.split(":", item)
        if len(arr1) >= 2:
            value = arr1[1]         # value = -0.007
            print(value)
            # 使用正则表达式进行匹配和提取
            pattern = r'\((.*?)\]'  # 匹配以括号开头，中间任意字符，以]结尾的部分
            match = re.search(pattern, arr1[0])
            if match:
                condition_text = match.group(1)  # 获取匹配结果中的第一个捕获组（codition_text: 0,11; 11,20; 20,32）
                # print(condition_text)
                arr3 = re.split(",", condition_text)
                print(arr3)
                if len(arr3)<2:
                    print("分段函数参数错误!");
                    continue
                if height_km > float(arr3[0]) and height_km <= float(arr3[1]):
                    return value
            else:
                print("No match found.")
                continue
        else:
            print("分段函数参数错误！")
            continue
    return 0


def testGetParamByHeight():
    param1 =  "(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001"
    height = "22000"
    ret = getParamByHeight(height, param1)
    print(ret)
def testParseParamStr():
    # param1 =  "(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001"
    param1 = "(0,50]km:3.5;(50,70]km:0.35"
    param2 = "(1,1000]Pa:10;(1000,110000]Pa:50"
    param3 = "(0,45]km:0.8;(45,55]km:1.0;(55,60]km:1.5;(60,70]km:2.3"
    param4 = "(0,30]km:0.1;(30,40]km:0.6;(40,45]km:1.3;(45,50]km:3.0;(50,55]km:5.1;(55,60]km:6.3;(60,70]km:9.0"
    param5 = "(0,50]km:0.1;(50,58]km:0.4;(58,60]km:0.9;(60,70]km:6.0"
    list = [param1, param2, param3, param4, param5]
    for item in list:
        print("=====")
        ret = parseParamStr(item)
        print(item)
        print(ret)


def parse_QQ_data(filepath):
    if not os.path.exists(filepath):
        print('File Error(not existed): {}'.format(filepath))
        return 0
    head_index = get_head_index_QQ(filepath)
    filename = os.path.basename(filepath)
    fileinfo = filename.split('_')
    device_code = fileinfo[1]       # 设备编码：探空气球等
    data_type = fileinfo[2]

    data = np.empty(0)
    if data_type == "ATDW":
        encode = get_file_encod(filepath)
        # data = pd.read_csv(filepath, skiprows=head_index, encoding=encode, sep='\s+').values
        probe_data = pd.read_csv(filepath, skiprows=head_index, encoding=encode, sep='\s+')
        data = probe_data.values
        # 原位火箭数据(s7)
        # 原位火箭错位head
        #      0           1          2               3                4              5              6                7                8                  9
        # ['Row_Number', 'Time', 'Longitude(°)', 'Latitude(°)', 'Altitude(m)', 'Atmospheric', 'Pressure(Pa)', 'Temperature(℃)', 'Density(kg/m³)', 'Velocity(m/s)', 'Wind_angle(°)']
        #     行号        时间      经度             纬度            高度            气压          温度          密度                  风速               风向
        #      1         14:00:00   87.380898      42.246109         11970         0.00254         97.35        1.195160                300                9.6            NaN
        data = data[:, (4, 5, 6, 7, 8, 9)]      # 0高度,1气压,2温度,3密度,4风速,5风向
    else:
        print('Filename Error: {}'.format(filepath))

    # plt.figure()
    # plt.plot(data[:, 1:], data[:, 0])
    # plt.ylabel('Height/km')
    # plt.show()
    data[data == -9999] = np.nan
    return data


def get_head_index_QQ(filepath):
    encode = get_file_encod(filepath)
    with open(filepath, 'r', encoding=encode) as f:
        filehead_all = f.readlines()
    head_index = 0
    while not (filehead_all[head_index].strip().startswith('Rang') or filehead_all[head_index].strip().startswith('Alt') or \
               filehead_all[head_index].strip().startswith('#End of Header') or filehead_all[head_index].strip().startswith('高度(距地)') or \
               filehead_all[head_index].strip().startswith('Row_Number')):
        head_index += 1
    return head_index


def get_file_encod(file_path):
    # 检查文件编码
    encode = get_file_encoding(file_path)
    print("file_encode:", encode)
    if encode == "utf-8":
        read_file_encoding = 'utf8'
    else:
        read_file_encoding = 'GBK'
    return read_file_encoding


if __name__ == "__main__":
    L2_filepath = r"C:\tmp\test_data_From_YJ_2311\原位探测new\探空火箭数据-1个文件\YL_BRS_ATDW_14DS20220101_20240125140000_V01.01.DAT"
    # L0_path = r"D:\1.Work\1.Data\新型地基\数据\实测数据\20230920\2023年09月20日数据分享\L0"
    ret = evaluate_uncertainty_rocket(L2_filepath, {})
    # ret = evaluate_uncertainty_rocket(L2_filepath, {}, 1.0, 'T')
    print("======================")
    print(ret)

    # 测试
    # testGetParamByHeight()
    # testParseParamStr()